clear all

% Simulate (gamma = 2 x estimated value), partial equilibrium

global Dividend

Dividend=0;

% Given Par2 optimal:
n = 50; tau = 0.1;
r = 0.036;
Par1 = zeros(31,1);
load ParEst_eh0_time0_Mex ParEst_eh0_time0_Mex
load ParEst_eh0_time1_Mex ParEst_eh0_time1_Mex
Par1(1:29) = ParEst_eh0_time0_Mex;
Par1(21) = 10^Par1(21);
Par1(22:27) = Par1(22:27)*1000;
Par1(30) = 2*500*ParEst_eh0_time1_Mex(1); 
Par1(31) = 2*500*ParEst_eh0_time1_Mex(2); 

load matrix_time1_age0_eh0_comp1.out
% Health transitions and unhealthy proportion
phl_data=matrix_time1_age0_eh0_comp1(1);
plh_data=matrix_time1_age0_eh0_comp1(2);
prob_l=matrix_time1_age0_eh0_comp1(3);

load W_time1_age0_eh0_comp11.out
w_f_g=W_time1_age0_eh0_comp11; % Read wages, g and f distributions
wf_1=w_f_g(:,1);        % Wage in formal sector
wi_1=w_f_g(:,2);        % Wage in informal sector
wf_2=w_f_g(:,7);      % Wage in formal sector
wi_2=w_f_g(:,8);      % Wage in informal sector

[Ff_1,ff_1,Fi_1,fi_1,Ff_2,ff_2,Fi_2,fi_2]=solve_wage_offer_dist(Par1,wf_1,wi_1,wf_2,wi_2,n);
[Wff_h,Wfi_h,Wfn_h,Wif_h,Wnf_h,Wii_h,Win_h,Wni_h,Wnn_h,Wff_l,Wfi_l,Wfn_l,Wif_l,Wnf_l,Wii_l,Win_l,Wni_l,Wnn_l,Rw_fi_ni_l,Rw_ii_ni_l,Rw_ii_in_l,Rw_if_in_l,...
    Rw_fi_ni_h,Rw_ii_ni_h,Rw_ii_in_h,Rw_if_in_h]=solve_value_function_A1_1(Par1,Ff_1,wf_1,Fi_1,wi_1,Ff_2,wf_2,Fi_2,wi_2,n);
[HHstates,wage_moments,transitions,phl,plh,gff,gfi,gif,gii,gfn,gnf,gin,gni]=data_simulation_A12_end(Par1,ff_1,wf_1,fi_1,wi_1,ff_2,wf_2,fi_2,wi_2,Wff_h,Wfi_h,Wfn_h,Wif_h,Wnf_h,Wii_h,Win_h,Wni_h,Wnn_h,Wff_l,Wfi_l,Wfn_l,Wif_l,Wnf_l,Wii_l,Win_l,Wni_l,Wnn_l);
wage_moments_gamma_1_PE = wage_moments;
HH_states_gamma_1_PE = HHstates;
save wage_moments_gamma_1_PE wage_moments_gamma_1_PE
save HH_states_gamma_1_PE HH_states_gamma_1_PE


mff=HHstates(1+9);
mfi=HHstates(2+9);
mfn=HHstates(3+9);
mif=HHstates(4+9);
mnf=HHstates(5+9);
mii=HHstates(6+9);
min_=HHstates(7+9);
mni=HHstates(8+9);
mnn=HHstates(9+9);
Un_Head = mnf + mni + mnn; Un_Spouse = mnn + min_ + mfn; Inf_rate = (mnn + min_ + mni + mii);
Agg_stocks_gamma_1_PE = [Inf_rate; Un_Head; Un_Spouse];
gf_1 = zeros (50,1); gf_2 = zeros (50,1); gi_1 =  zeros (50,1); gi_2 = zeros (50,1);
 for i = 1:50
 gf_1(i) =  (mff*(sum(gff(i,:))) + mfi*(sum(gfi(i,:))) + mfn*gfn(i))/(mff+mfi+mfn);   
 end

  for i = 1:50
 gi_1(i) =  (mif*(sum(gif(i,:))) + mii*(sum(gii(i,:))) + min_*gin(i))/(mif+mii+min_);   
  end
 
   for i = 1:50
 gf_2(i) =  (mff*(sum(gff(:,i))) + mif*(sum(gif(:,i))) + mnf*gnf(i))/(mff+mif+mnf);   
   end
 
    for i = 1:50
 gi_2(i) =  (mfi*(sum(gfi(:,i))) + mii*(sum(gii(:,i))) + mni*gni(i))/(mfi+mii + mni);   
 end
save Agg_stocks_gamma_1_PE Agg_stocks_gamma_1_PE

% have to construct the joint probabilities (for ex, gff(wf_1,wf_2)) or integrate over the
% marginal distributions gf_1 and gf_2, whichever is easier..

mff_h=HHstates(1);
mfi_h=HHstates(2);
mfn_h=HHstates(3);
mif_h=HHstates(4);
mnf_h=HHstates(5);
mii_h=HHstates(6);
min_h=HHstates(7);
mni_h=HHstates(8);
mnn_h=HHstates(9);
mff_l=(mff- mff_h*(1-prob_l))/ prob_l;
mfi_l=(mfi- mfi_h*(1-prob_l))/ prob_l;
mfn_l=(mfn- mfn_h*(1-prob_l))/ prob_l;
mif_l=(mif- mif_h*(1-prob_l))/ prob_l;
mnf_l=(mnf- mnf_h*(1-prob_l))/ prob_l;
mii_l=(mii- mii_h*(1-prob_l))/ prob_l;
min_l=(min_- min_h*(1-prob_l))/ prob_l;
mni_l=(mni- mni_h*(1-prob_l))/ prob_l;
mnn_l=(mnn- mnn_h*(1-prob_l))/ prob_l;

Welfare_total_h=r*(sum(sum(mff_h*Wff_h.*gff+mfi_h*Wfi_h.*gfi+mif_h*Wif_h.*gif+mii_h*Wii_h.*gii)))+r*sum(mfn_h*Wfn_h.*gfn+mnf_h*Wnf_h.*gnf+min_h*Win_h.*gin+mni_h*Wni_h.*gni)+r*mnn_h*Wnn_h;
Welfare_total_l=r*(sum(sum(mff_l*Wff_l.*gff+mfi_l*Wfi_l.*gfi+mif_l*Wif_l.*gif+mii_l*Wii_l.*gii)))+r*sum(mfn_l*Wfn_l.*gfn+mnf_l*Wnf_l.*gnf+min_l*Win_l.*gin+mni_l*Wni_l.*gni)+r*mnn_l*Wnn_l;
Welfare_total=prob_l*Welfare_total_l+(1-prob_l)*Welfare_total_h;
Welfare_formal1_h=(r*(sum(sum(mff_h*Wff_h.*gff+mfi_h*Wfi_h.*gfi)))+r*sum(mfn_h*Wfn_h.*gfn))/(mfn_h+mfi_h+mff_h);
Welfare_informal1_h=(r*(sum(sum(mif_h*Wfi_h.*gif+mii_h*Wii_h.*gii)))+r*sum(min_h*Win_h.*gin))/(mif_h+mii_h+min_h);
Welfare_nonemployed1_h=(r*sum(mnf_h*Wnf_h.*gnf+mni_h*Wni_h.*gni)+r*mnn_h*Wnn_h)/(mnf_h+mni_h+mnn_h);
Welfare_formal2_h=(r*(sum(sum(mff_h*Wff_h.*gff+mif_h*Wif_h.*gif)))+r*sum(mnf_h*Wnf_h.*gnf))/(mnf_h+mif_h+mff_h);
Welfare_informal2_h=(r*(sum(sum(mfi_h*Wfi_h.*gfi+mii_h*Wii_h.*gii)))+r*sum(mni_h*Wni_h.*gni))/(mni_h+mii_h+mfi_h);
Welfare_nonemployed2_h=(r*sum(mfn*Wfn_h.*gfn+min_h*Win_h.*gin)+r*mnn_h*Wnn_h)/(min_h+mfn+mnn_h);
Welfare_formal1_l=(r*(sum(sum(mff_l*Wff_l.*gff+mfi_l*Wfi_l.*gfi)))+r*sum(mfn_l*Wfn_l.*gfn))/(mfn_l+mfi_l+mff_l);
Welfare_informal1_l=(r*(sum(sum(mif_l*Wfi_l.*gif+mii_l*Wii_l.*gii)))+r*sum(min_l*Win_l.*gin))/(mif_l+mii_l+min_l);
Welfare_nonemployed1_l=(r*sum(mnf_l*Wnf_l.*gnf+mni_l*Wni_l.*gni)+r*mnn_l*Wnn_l)/(mnf_l+mni_l+mnn_l);
Welfare_formal2_l=(r*(sum(sum(mff_l*Wff_l.*gff+mif_l*Wif_l.*gif)))+r*sum(mnf_l*Wnf_l.*gnf))/(mnf_l+mif_l+mff_l);
Welfare_informal2_l=(r*(sum(sum(mfi_l*Wfi_l.*gfi+mii_l*Wii_l.*gii)))+r*sum(mni_l*Wni_l.*gni))/(mni_l+mii_l+mfi_l);
Welfare_nonemployed2_l=(r*sum(mfn*Wfn_l.*gfn+min_l*Win_l.*gin)+r*mnn_l*Wnn_l)/(min_l+mfn+mnn_l);
Welfare_formal1=prob_l*Welfare_formal1_l+(1-prob_l)*Welfare_formal1_h;
Welfare_informal1=prob_l*Welfare_informal1_l+(1-prob_l)*Welfare_informal1_h;
Welfare_nonemployed1=prob_l*Welfare_nonemployed1_l+(1-prob_l)*Welfare_nonemployed1_h;
Welfare_formal2=prob_l*Welfare_formal2_l+(1-prob_l)*Welfare_formal2_h;
Welfare_informal2=prob_l*Welfare_informal2_l+(1-prob_l)*Welfare_informal2_h;
Welfare_nonemployed2=prob_l*Welfare_nonemployed2_l+(1-prob_l)*Welfare_nonemployed2_h;
Welfare_gamma_1_PE=[Welfare_total;Welfare_total_h;Welfare_total_l;Welfare_formal1;Welfare_informal1;Welfare_nonemployed1;Welfare_formal2;Welfare_informal2;Welfare_nonemployed2]; % IN THOUSANDS, IF NORMALIZED
save Welfare_gamma_1_PE Welfare_gamma_1_PE

% Generate LOG WAGES by the F distribution 
[wvector_f_1_est,Fvector_f_1_w]=main_quantiles(1-Ff_1,wf_1,n);
[wvector_i_1_est,Fvector_i_1_w]=main_quantiles(1-Fi_1,wi_1,n);
[wvector_f_2_est,Fvector_f_2_w]=main_quantiles(1-Ff_2,wf_2,n);
[wvector_i_2_est,Fvector_i_2_w]=main_quantiles(1-Fi_2,wi_2,n);
mean_wf_1f=log(sum(wf_1.*ff_1));
mean_wi_1f=log(sum(wi_1.*fi_1));
mean_wf_2f=log(sum(wf_2.*ff_2));
mean_wi_2f=log(sum(wi_2.*fi_2));
wage_moments_under_F=[[wvector_f_1_est;mean_wf_1f],[wvector_i_1_est;mean_wi_1f],[wvector_f_2_est;mean_wf_2f],[wvector_i_2_est;mean_wi_2f]];
save wage_moments_under_f_gamma_1_PE wage_moments_under_F

% Transitions 
Trans_head = [transitions(3); transitions(4); transitions(1); transitions(5); transitions(2); transitions(6); transitions(13); transitions(14)];
Trans_spouse = [transitions(9); transitions(10); transitions(7); transitions(11); transitions(8); transitions(12); transitions(15); transitions(16)];
Transitions_gamma_1_PE = [Trans_head,Trans_spouse];
save Transitions_gamma_1_PE Transitions_gamma_1_PE

% Constructing produtivities for simulations

% marginal distributions
mf_1=mff+mfi+mfn;
mi_1=mif+mii+min_;
mf_2=mff+mif+mnf;
mi_2=mfi+mii+mni;
% gf_1=gfn*mfn/mf_1+sum(gfi')'*mfi/mf_1+sum(gff')'*mff/mf_1;
% gi_1=gin*min_/mi_1+sum(gii')'*mii/mi_1+sum(gif')'*mif/mi_1;
% gf_2=gnf*mnf/mf_2+sum(gif)'*mif/mf_2+sum(gff)'*mff/mf_2;
% gi_2=gni*mni/mi_2+sum(gii)'*mii/mi_2+sum(gfi)'*mfi/mi_2;


% firm size
el_f_1=mf_1*gf_1./ff_1; 
% el_f_1(1)=el_f_1(2); % may need smoothing
% el_f_1(n)=el_f_1(n-1);
el_i_1=mi_1*gi_1./fi_1;
% el_i_1(1)=el_i_1(2); % may need smoothing
% el_i_1(n)=el_i_1(n-1);
el_f_2=mf_2*gf_2./ff_2; 
% el_f_2(1)=el_f_2(2); % may need smoothing
% el_f_2(n)=el_f_2(n-1);
el_i_2=mi_2*gi_2./fi_2;
% el_i_2(1)=el_i_2(2); % may need smoothing
% el_i_2(n)=el_i_2(n-1);

load fitted_model p_f_1 p_i_1 p_f_2 p_i_2 Ff_1 Fi_1 Ff_2 Fi_2 ff_1 fi_1 ff_2 fi_2

% M (total employed workforce) in 2005 = 41441076
% Fraction of firms in the formal sector in 2002 includes self-emp. (ENAMIN/INEGI: 23.7%)
% Number of formal firms (IMSS registry, BOSCH AND CAMPOS-VASQUEZ 2005): 800,000
% Number of firms = 800000/0.237
M_N=41441076/(800000/0.237);

pi_f_1=(p_f_1-wf_1*(1+tau)).*el_f_1*M_N/0.237;
pi_i_1=(p_i_1-wi_1).*el_i_1*M_N/0.763;
pi_f_2=(p_f_2-wf_2*(1+tau)).*el_f_2*M_N/0.237;
pi_i_2=(p_i_2-wi_2).*el_i_2*M_N/0.763;
% may need smoothing
% pi_f_1s=smoothdata(pi_f_1(1:25),'rlowess',20);
% pi_i_1s=smoothdata(pi_i_1(1:18),'rlowess',20);
% pi_f_2s=smoothdata(pi_f_2(1:25),'rlowess',20);
% pi_i_2s=smoothdata(pi_i_2(1:18),'rlowess',20);

profitFormalFirms_1=1./(M_N*mf_1/0.237)*sum(pi_f_1(1:25).*ff_1(1:25)); % per worker 
profitInformalFirms_1=1./(M_N*mi_1/0.763)*sum(pi_i_1(1:18).*fi_1(1:18)); % per worker 
totalprofitFirms_1=mf_1*profitFormalFirms_1+mi_1*profitInformalFirms_1;  
profitFormalFirms_2=1./(M_N*mf_2/0.237)*sum(pi_f_2(1:25).*ff_2(1:25)); % per worker 
profitInformalFirms_2=1./(M_N*mi_2/0.763)*sum(pi_i_2(1:18).*fi_2(1:18)); % per worker 
totalprofitFirms_2=mf_2*profitFormalFirms_2+mi_2*profitInformalFirms_2;  
Profit_gamma_1_PE=[totalprofitFirms_1;profitFormalFirms_1;profitInformalFirms_1;NaN;totalprofitFirms_2;profitFormalFirms_2;profitInformalFirms_2]; % per worker
save Profit_gamma_1_PE Profit_gamma_1_PE

% normalized firm size
[elvector_f_1,Fvector_f_1]=main_quantiles2(1-Ff_1,el_f_1,n)
[elvector_i_1,Fvector_i_1]=main_quantiles2(1-Fi_1,el_i_1,n)
[elvector_f_2,Fvector_f_2]=main_quantiles2(1-Ff_2,el_f_2,n)
[elvector_i_2,Fvector_i_2]=main_quantiles2(1-Fi_2,el_i_2,n)
mean_el_f_1=sum(el_f_1.*ff_1);
mean_el_i_1=sum(el_i_1.*fi_1);
mean_el_f_2=sum(el_f_2.*ff_2);
mean_el_i_2=sum(el_i_2.*fi_2);
Firm_size_gamma_1_PE=[[elvector_f_1,elvector_i_1,elvector_f_2,elvector_i_2];[mean_el_f_1,mean_el_i_1,mean_el_f_2,mean_el_i_2]];
save Firm_size_gamma_1_PE Firm_size_gamma_1_PE

